Takehome_Ex01

1 Setting the Scene

Road traffic accidents are a significant global issue, causing approximately 1.19 million deaths and leaving 20 to 50 million people with non-fatal injuries annually, according to the World Health Organization (WHO). Vulnerable road users, like pedestrians, cyclists, and motorcyclists, make up more than half of these fatalities, with most incidents occurring in low- and middle-income countries. Thailand, in particular, has one of the highest road traffic death rates in Southeast Asia, with around 20,000 deaths each year. Economic impacts are also substantial, costing countries about 3% of their GDP. In Thailand, national highways see the highest concentration of accidents, especially in straight road sections and other accident-prone zones like intersections and curves.

2 Objectives

3 Data Collection

3.1 Loading Packages

pacman::p_load(sf, tidyverse, tmap, spNetwork, sp,spatstat ,dplyr)

3.2 Data Collect

Three data sets are used in this exercise, they are:

1. Thailand Road Accident [2019-2022] on Kaggle. This dataset offers comprehensive statistics on road accidents recorded in Thailand from approximately 2019 to 2022, covering various aspects of the incidents.

2. Thailand Roads (OpenStreetMap Export) on HDX. This dataset, sourced from OpenStreetMap, provides a detailed map of Thailand’s road network.

3. Thailand - Subnational Administrative Boundaries on HDX. This dataset provides comprehensive geographic data on Thailand’s subnational administrative boundaries, covering provinces, districts, and subdistricts.

4 KDE analysis

In this part, we will conduct a kernel density analysis for the Bangkok Metropolitan Region (BMR) to simply visualize and identify areas with a higher concentration of car accidents.

4.1 Data Preparation

4.4.1 Data import

acc <- read.csv("data/aspatical/thai_road_accident_2019_2022.csv") 

tr <- st_read(dsn="data/geospatial", layer = 'hotosm_tha_roads_lines_shp')

Reading layer `hotosm_tha_roads_lines_shp' from data source 
  `D:\FuWanqian\ISSS608-VAA\Takehome_Ex\Takehome_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 2792361 features and 14 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 97.34457 ymin: 5.643645 xmax: 105.6528 ymax: 20.47168
Geodetic CRS:  WGS 84

ad <- st_read(dsn="data/geospatial", layer = 'tha_admbnda_adm1_rtsd_20220121')

Reading layer `tha_admbnda_adm1_rtsd_20220121' from data source 
  `D:\FuWanqian\ISSS608-VAA\Takehome_Ex\Takehome_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 77 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 97.34336 ymin: 5.613038 xmax: 105.637 ymax: 20.46507
Geodetic CRS:  WGS 84

Here we first import the administrative level 1 data, which includes the provincial boundaries of Thailand.

4.4.2 Data selection

bmr_acc<- acc %>%
  filter(province_en %in% c("Bangkok", "Nakhon Pathom", "Pathum Thani", 
                        "Nonthaburi", "Samut Prakan", "Samut Sakhon"))

bmr_ad<- ad %>%
  filter(ADM1_EN %in% c("Bangkok", "Nakhon Pathom", "Pathum Thani", 
                        "Nonthaburi", "Samut Prakan", "Samut Sakhon"))

Here we filter the BMR data from both the accident data and the provincial boundary data of Thailand.”

4.4.3 Drop missing value for accident data

missing_coords <- bmr_acc[is.na(bmr_acc$longitude) | is.na(bmr_acc$latitude), ]

missing_coords_count <- nrow(missing_coords)

missing_coords_count

[1] 350

The number of missing coordinates is 350, so using below codes to drop missing value.

bmr_acc2 <- bmr_acc[!is.na(bmr_acc$longitude) & !is.na(bmr_acc$latitude), ]

4.4.4 Transform data format

It is important for us to ensure that all the data are projected in same projection system, then here we transform BMR accident data and BMR boundary to EPSG:32647 (UTM Zone 47N).

bmr_acc_data <- st_as_sf(bmr_acc2, coords = c("longitude", "latitude"), crs = 4326)

bmr_acc_data_utm <- st_transform(bmr_acc_data, crs = 32647)
bmr_acc_data_utm

Simple feature collection with 12986 features and 16 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 591277.5 ymin: 1486846 xmax: 710166.1 ymax: 1576520
Projected CRS: WGS 84 / UTM zone 47N
First 10 features:
   acc_code incident_datetime report_datetime  province_th   province_en
1    571882     2019/1/1 2:25  2019/1/2 17:32       นครปฐม Nakhon Pathom
2    600001     2019/1/1 3:00  2019/1/5 10:33        นนทบุรี    Nonthaburi
3    605043     2019/1/1 3:00  2019/3/29 8:22   สมุทรปราการ  Samut Prakan
4    629691     2019/1/1 3:05   2019/1/1 3:05 กรุงเทพมหานคร       Bangkok
5    571887     2019/1/1 4:30  2019/1/2 17:32       นครปฐม Nakhon Pathom
6    599234     2019/1/1 4:45   2019/1/2 8:28   สมุทรปราการ  Samut Prakan
7    599990     2019/1/1 5:30  2019/1/5 10:35     สมุทรสาคร  Samut Sakhon
8    612045     2019/1/1 5:30 2019/10/9 13:50        นนทบุรี    Nonthaburi
9    629689     2019/1/1 5:42   2019/1/1 5:42 กรุงเทพมหานคร       Bangkok
10   607046     2019/1/1 6:30 2019/5/15 13:24       ปทุมธานี  Pathum Thani
                             agency
1         department of rural roads
2            department of highways
3            department of highways
4  expressway authority of thailand
5         department of rural roads
6            department of highways
7            department of highways
8            department of highways
9  expressway authority of thailand
10           department of highways
                                            route          vehicle_type
1   แยกทางหลวงหมายเลข 4 (กม.ที่ 51+920) - บ้านวัดละมุด            motorcycle
2                            คลองวัดแดง - บางบัวทอง private/passenger car
3                         ราษฎร์บูรณะ - พระสมุทรเจดีย์ private/passenger car
4                                    บางพลี-สุขสวัสดิ์                 other
5  แยกทางหลวงหมายเลข 3310 (กม.ที่ 10+700) - บ้านกลาง            motorcycle
6                                   บางพลี - บางบ่อ            motorcycle
7                แสมดำ - สะพานข้ามแม่น้ำท่าจีนฝั่งตะวันตก            motorcycle
8                     แยกพงษ์เพชร - สะพานพระนั่งเกล้า private/passenger car
9                                      บางนา-ชลบุรี                 other
10                               ปทุมธานี - ท้ายเกาะ private/passenger car
                                  presumed_cause
1                                       speeding
2                                       speeding
3             running red lights/traffic signals
4                                          other
5                                       speeding
6         driving under the influence of alcohol
7                                       speeding
8  cutting in closely by people/vehicles/animals
9                                          other
10                                      speeding
                      accident_type number_of_vehicles_involved
1  rollover/fallen on straight road                           1
2  rollover/fallen on straight road                           1
3  collision at intersection corner                           2
4                             other                           1
5  rollover/fallen on straight road                           1
6  rollover/fallen on straight road                           1
7  rollover/fallen on straight road                           1
8  rollover/fallen on straight road                           1
9                             other                           1
10 rollover/fallen on straight road                           1
   number_of_fatalities number_of_injuries weather_condition road_description
1                     0                  2             clear    straight road
2                     0                  1             clear    straight road
3                     0                  0             clear            other
4                     0                  1             clear            other
5                     0                  1             clear    straight road
6                     1                  0             clear    straight road
7                     1                  0             clear    straight road
8                     0                  0             clear    straight road
9                     0                  0             clear    straight road
10                    0                  1             clear    straight road
   slope_description                 geometry
1           no slope POINT (627012.3 1533381)
2           no slope POINT (655411.7 1531940)
3              other POINT (670903.3 1503427)
4              other POINT (674117.8 1525144)
5           no slope POINT (641347.8 1526237)
6           no slope POINT (693488.9 1502876)
7           no slope POINT (640469.3 1501125)
8           no slope POINT (661585.2 1533339)
9              other POINT (674117.8 1525144)
10          no slope POINT (664551.7 1554801)

Simple feature collection with 12986 features and 16 fields
Geometry type: POINT Dimension: XY
Bounding box: xmin: 591277.5 ymin: 1486846 xmax: 710166.1 ymax: 1576520
Projected CRS: WGS 84 / UTM zone 47N

bmr_ad_data_utm <- st_transform(bmr_ad, crs = 32647)
bmr_ad_data_utm

Simple feature collection with 6 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 587893.5 ymin: 1484414 xmax: 712440.5 ymax: 1579076
Projected CRS: WGS 84 / UTM zone 47N
  Shape_Leng Shape_Area       ADM1_EN      ADM1_TH ADM1_PCODE ADM1_REF
1   2.417227 0.13133873       Bangkok กรุงเทพมหานคร       TH10     <NA>
2   1.695100 0.07926199  Samut Prakan   สมุทรปราการ       TH11     <NA>
3   1.251111 0.05323766    Nonthaburi        นนทบุรี       TH12     <NA>
4   1.884945 0.12698345  Pathum Thani       ปทุมธานี       TH13     <NA>
5   2.463030 0.17891420 Nakhon Pathom       นครปฐม       TH73     <NA>
6   1.566369 0.07155983  Samut Sakhon     สมุทรสาคร       TH74     <NA>
  ADM1ALT1EN ADM1ALT2EN ADM1ALT1TH ADM1ALT2TH  ADM0_EN   ADM0_TH ADM0_PCODE
1       <NA>       <NA>       <NA>       <NA> Thailand ประเทศไทย         TH
2       <NA>       <NA>       <NA>       <NA> Thailand ประเทศไทย         TH
3       <NA>       <NA>       <NA>       <NA> Thailand ประเทศไทย         TH
4       <NA>       <NA>       <NA>       <NA> Thailand ประเทศไทย         TH
5       <NA>       <NA>       <NA>       <NA> Thailand ประเทศไทย         TH
6       <NA>       <NA>       <NA>       <NA> Thailand ประเทศไทย         TH
        date    validOn    validTo                       geometry
1 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((674339.8 15...
2 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((687139.8 15...
3 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((644817.9 15...
4 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((704086 1575...
5 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((631987.6 15...
6 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((641549.1 15...

Simple feature collection with 6 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 587893.5 ymin: 1484414 xmax: 712440.5 ymax: 1579076
Projected CRS: WGS 84 / UTM zone 47N

4.2 Mapping the data

plot(st_geometry(bmr_ad_data_utm), col = "lightgrey", border = "black", main = "BMR Accident Locations (Base R)")

plot(st_geometry(bmr_acc_data_utm), col = "black", pch = 19, cex = 0.1, add = TRUE)  

tm_shape(bmr_ad_data_utm) +
  tm_borders(col = "black", lwd = 1) +   borders

  tm_shape(bmr_acc_data_utm) +
  tm_dots(col = "black", size = 0.01) +  
  tm_layout(title = "BMR Accident Locations (tmap)")

The first plot shows the locations of traffic accidents in the Bangkok Metropolitan Region, with black dots concentrated in the city center and along major roads, especially near intersections and busy traffic areas. Accident density varies across different provinces, with hotspots mainly distributed along major roads.

4.3 Convert data to ppp format

4.3.1 sp format

bmr_acc_data_sp <- as_Spatial(bmr_acc_data_utm)
bmr_boundary_sp <- as_Spatial(bmr_ad_data_utm)

4.3.2 ppp format

acc_coords <- coordinates(bmr_acc_data_sp)
bbox_values <- bbox(bmr_boundary_sp)

bmr_window <- owin(xrange = c(bbox_values[1, 1], bbox_values[1, 2]), 
                   yrange = c(bbox_values[2, 1], bbox_values[2, 2]))

bmr_acc_ppp <- ppp(x = acc_coords[, 1], y = acc_coords[, 2], window = bmr_window)
plot(bmr_acc_ppp, main = "BMR Accident Locations", cex = 0.5)

4.4 Handling duplicate points

4.4.1 check duplicates

sum(multiplicity(bmr_acc_ppp) > 1)

[1] 2293

4.4.2 Apply jitter to handle duplicates

bmr_acc_ppp_jit <- rjitter(bmr_acc_ppp, retry = TRUE, nsim = 1, drop = TRUE)
any(duplicated(bmr_acc_ppp_jit))

[1] FALSE

4.5 Combine point events object and owin object

bmr_acc_ppp_final <- bmr_acc_ppp_jit[bmr_window]

4.6 Kernel Density Estimation (KDE) Analysis

sigma_value <- 12000  
kde_bmr_acc <- density(bmr_acc_ppp_final,  
                       sigma = sigma_value, 
                       edge = TRUE,  
                       kernel = "gaussian")  

plot(kde_bmr_acc, main = "KDE of Accident Data")
plot(st_geometry(bmr_ad_data_utm), add = TRUE, border = "black", lwd = 2)

From the map, it is clear that central Bangkok and parts of Samut Prakan exhibit the highest density of car accidents, indicating that accidents are more concentrated in these urban areas. In contrast, the outer regions such as Nakhon Pathom, Pathum Thani, and Nonthaburi show much lower accident densities.

4.7 G Function

Null Hypothesis (H0): The spatial distribution of car accidents follows a completely spatially random (CSR) pattern, meaning that accidents occur independently and uniformly over the study area.

Alternative Hypothesis (H1): The spatial distribution of car accidents is not random and exhibits clustering or dispersion, implying that accidents are more likely to occur near each other or further apart than expected under CSR.

G_CK = Gest(bmr_acc_ppp_final, correction = "border")
plot(G_CK, xlim=c(0,500))

The plot shows that the observed G function (black line) rises more quickly than the expected Poisson function (red dashed line), indicating that car accidents are spatially clustered.

5 NKDE analysis

Based on the KDE analysis, we found that Bangkok has the highest concentration of car accidents. Therefore, we refine the NKDE (Network Kernel Density Estimation) analysis to focus specifically on Bangkok’s road network, allowing for a more precise identification of accident hotspots along key roads and intersections.

5.1 Data select

Lots of side roads in the road data may influence the accuracy of the NKDE results, so here we filter the dataset to include only major roads.

main_rd <- tr %>%
  filter(highway %in% c("motorway", "trunk", "primary", "secondary"))
saveRDS(main_rd,"data/geospatial/main_rd.rds")

Here we select bangkok city boundary.

selected_cities <- c("Bangkok")

bangkok_city <- bmr_ad_data_utm %>%
  filter(ADM1_EN %in% selected_cities)

saveRDS(bangkok_city,"data/geospatial/bmr_city.rds")

Then we select the main roads within Bangkok by intersecting the road and city boundary datasets.

main_rd_format <- st_transform(main_rd, crs = 32647)
saveRDS(main_rd_format,"data/geospatial/main_rd_format.rds")

main_bangkok_roads <- st_intersection(main_rd_format, bangkok_city)

saveRDS(main_bangkok_roads, "data/geospatial/main_bangkok_roads2.rds")

main_bangkok_roads <- readRDS("data/geospatial/main_bangkok_roads2.rds")

Filters the accident data to include only records from the Bangkok province.

bangkok_acc <- bmr_acc_data_utm %>%
  filter(province_en %in% "Bangkok")

5.2 Transfer to linestring

main_bangkok_roads <- main_bangkok_roads %>%
  filter(st_geometry_type(main_bangkok_roads) %in% c("LINESTRING", "MULTILINESTRING"))


main_bangkok_roads <- st_cast(main_bangkok_roads, "LINESTRING", group_or_split = TRUE)

5.3 Generate lixels

lixels_bangkok <- lixelize_lines(main_bangkok_roads,
                         10000,        
                         mindist = 5000) 


samples_bangkok <- lines_center(lixels_bangkok)

5.4 Perform NKDE

nkde_result_bangkok <- nkde(
  lines = lixels_bangkok,                      
  events = bangkok_acc,                     
  w = rep(1, nrow(bangkok_acc)),            
  samples = samples_bangkok,                   
  kernel_name = "quartic",                     
  bw = 500,                                    
  div = "bw",                                  
  method = "simple",                          
  grid_shape = c(200, 200),                    
  verbose = TRUE                               
)

nkde_result_bangkok <- readRDS("data/aspatical/nkde_result_bangkok.rds")

head(nkde_result_bangkok, 10)

 [1] 0.000000e+00 0.000000e+00 0.000000e+00 2.435628e-06 0.000000e+00
 [6] 0.000000e+00 3.205377e-06 0.000000e+00 2.432477e-06 0.000000e+00

5.5 Visualize NKDE results

samples_bangkok$density <- nkde_result_bangkok
lixels_bangkok$density <- nkde_result_bangkok

samples_bangkok$density <- samples_bangkok$density * 10000
lixels_bangkok$density <- lixels_bangkok$density * 10000

tmap_mode('view')

tm_shape(lixels_bangkok) +
  tm_lines(col = "density", palette = "YlOrRd", title.col = "Density", lwd = 2) +
  tm_shape(bangkok_acc) +
  tm_dots(size = 0.1, col = "blue", alpha = 0.5, title = "Accidents")

tmap_mode('plot')

From above plot, the traffic accidents are almost entirely concentrated along the road network, particularly in the main roads and intersections of central Bangkok. The highest accident density is found along the main roads in central Bangkok, particularly in areas like Phra Nakhon, Pathum Wan, and Sathon. These indicate high-risk zones for accidents. In contrast, outer areas such as Nonthaburi and Samut Prakan have much lower accident densities, suggesting fewer incidents, likely due to lower traffic volumes or more dispersed road networks.

6 Spatio-Temporal analysis

In this part, we will analyze the accident distribution of BMR from a time dimension in order to understand how the frequency and location of accidents vary over time, helping to identify patterns and trends that could inform better traffic management and accident prevention strategies.

6.1 Split time

Here we divide the day into six time intervals.

bmr_acc$incident_datetime <- as.POSIXct(bmr_acc$incident_datetime, format="%Y/%m/%d %H:%M")

bmr_acc$hour <- hour(bmr_acc$incident_datetime)

bmr_acc$time_period <- cut(bmr_acc$hour, 
                           breaks = seq(0, 24, by = 2), 
                           include.lowest = TRUE, 
                           labels = FALSE)

accident_count <- bmr_acc %>%
  group_by(time_period) %>%
  summarise(count = n())

time_labels <- c("00:00-02:00", "02:00-04:00", "04:00-06:00", "06:00-08:00", 
                 "08:00-10:00", "10:00-12:00", "12:00-14:00", "14:00-16:00", 
                 "16:00-18:00", "18:00-20:00", "20:00-22:00", "22:00-24:00")


ggplot(accident_count, aes(x = factor(time_period, labels = time_labels), y = count)) +
  geom_line(group=1, color="blue") +
  geom_point(color="red") +
  labs(title = "Car Accident Count by Time of Day", x = "Time Period", y = "Accident Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

From the line plot above, we can see that the 08:00-10:00 has the highest accident count and 22:00-24:00 has the lowest count. Therefore, we will focus on these two time intervals for further analysis.

6.2 Time KDE analysis for BMR

6.2.1 create new data frame

bmr_acc$hour <- as.numeric(format(bmr_acc$incident_datetime, "%H"))


bmr_acc$time_period_flag <- ifelse(bmr_acc$hour >= 8 & bmr_acc$hour < 10, "08:00-10:00", 
                            ifelse(bmr_acc$hour >= 22 & bmr_acc$hour < 24, "22:00-24:00", NA))


filtered_acc <- bmr_acc %>%
  filter(!is.na(time_period_flag))

peak_acc <- filtered_acc %>%
  filter(time_period_flag == "08:00-10:00")

non_peak_acc <- filtered_acc %>%
  filter(time_period_flag == "22:00-24:00")

6.2.2 Transform accident data to EPSG:32647 (UTM Zone 47N)

peak_acc_clean <- peak_acc %>%
  filter(!is.na(longitude) & !is.na(latitude))

peak_acc_sf <- st_as_sf(peak_acc_clean, coords = c("longitude", "latitude"), crs = 4326)

peak_acc_utm <- st_transform(peak_acc_sf, crs = 32647)

non_peak_acc_clean <- non_peak_acc %>%
  filter(!is.na(longitude) & !is.na(latitude))

non_peak_acc_sf <- st_as_sf(non_peak_acc_clean, coords = c("longitude", "latitude"), crs = 4326)

non_peak_acc_utm <- st_transform(non_peak_acc_sf, crs = 32647)

6.2.3 Convert data to ppp

peak_acc_utm_sp <- as_Spatial(peak_acc_utm)
nonpeak_acc_utm_sp <- as_Spatial(non_peak_acc_utm)

bmr_boundary_sp <- as_Spatial(bmr_ad_data_utm)

peak_acc_coords <- coordinates(peak_acc_utm_sp)
nonpeak_acc_coords <- coordinates(nonpeak_acc_utm_sp)

bbox_values <- bbox(bmr_boundary_sp)

bmr_window <- owin(xrange = c(bbox_values[1, 1], bbox_values[1, 2]), 
                   yrange = c(bbox_values[2, 1], bbox_values[2, 2]))

peak_bmr_acc_ppp <- ppp(x = peak_acc_coords[, 1], y = peak_acc_coords[, 2], window = bmr_window)

nonpeak_bmr_acc_ppp <- ppp(x = nonpeak_acc_coords[, 1], y = nonpeak_acc_coords[, 2], window = bmr_window)

6.2.4 Handling duplicate points

peak_bmr_acc_jit <- rjitter(peak_bmr_acc_ppp, retry = TRUE, nsim = 1, drop = TRUE)

nonpeak_bmr_acc_jit <- rjitter(nonpeak_bmr_acc_ppp, retry = TRUE, nsim = 1, drop = TRUE)

Check for duplicates after jittering

any(duplicated(peak_bmr_acc_jit))

[1] FALSE

any(duplicated(nonpeak_bmr_acc_jit))

[1] FALSE

6.2.5 Combine point events object and owin object

peak_bmr_acc_final <- peak_bmr_acc_jit[bmr_window]
nonpeak_bmr_acc_final <- nonpeak_bmr_acc_jit[bmr_window]

6.2.6 Kernel Density Estimation (KDE) Analysis

sigma_value <- 5000  


peak_kde_bmr_acc <- density(peak_bmr_acc_final,  
                       sigma = sigma_value,  
                       edge = TRUE,  
                       kernel = "gaussian")  

sigma_value <- 5000  



nonpeak_kde_bmr_acc <- density(nonpeak_bmr_acc_final,  #
                       sigma = sigma_value,  
                       edge = TRUE,  
                       kernel = "gaussian")  

plot(peak_kde_bmr_acc, main = "KDE of Accident Data(peak)")
plot(st_geometry(bmr_ad_data_utm), add = TRUE, border = "black", lwd = 2)

plot(nonpeak_kde_bmr_acc, main = "KDE of Accident Data(non-peak)")
plot(st_geometry(bmr_ad_data_utm), add = TRUE, border = "black", lwd = 2)

In the peak image, we can see the accident hotspots are mainly concentrated in the city center of Bangkok and surrounding areas. These areas likely experience a higher frequency of accidents due to the heavy traffic flow and commuter pressure. Additionally, in the western part of Bangkok, accident density is also significant. This indicates that during peak hours, accidents are not only concentrated in the core city but also extend to outer areas connected to the city, where traffic pressure is also high.

In the non-peak period image, the distribution of accidents appears more dispersed, though the central and surrounding areas of Bangkok still remain hotspots. Compared to the peak period, the frequency of accidents has reduced, and the areas of concentration have become smaller. In non-peak times, the western part of the city sees a decrease in accident density, with fewer hotspots overall, particularly in areas farther from the core city. The smaller yellow regions indicate fewer accidents when the traffic flow is lower, showing that during non-peak hours, although accidents still occur, they are less frequent, and their distribution is more spread out across the region.

6.3 Time KDE for bangkok

From the time KDE for BMR, we can find that Bangkok has a higher car accident density in both time intervals. Therefore, we narrow our focus to Bangkok to examine the city’s accident distribution in more detail.

6.3.1 select data

acc data

peak_bangkok_acc <- peak_acc_utm %>%
  filter(province_en %in% "Bangkok")

nonpeak_bangkok_acc <- non_peak_acc_utm %>%
  filter(province_en %in% "Bangkok")

bangkok district data

bangkok_district <- st_read(dsn = "data/geospatial",
                        layer = "tha_admbnda_adm2_rtsd_20220121")%>%
  filter(ADM1_EN == "Bangkok")

Reading layer `tha_admbnda_adm2_rtsd_20220121' from data source 
  `D:\FuWanqian\ISSS608-VAA\Takehome_Ex\Takehome_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 928 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 97.34336 ymin: 5.613038 xmax: 105.637 ymax: 20.46507
Geodetic CRS:  WGS 84

bangkok_district_format <- st_transform(bangkok_district, crs = 32647)

6.3.2 Convert data to ppp

peak_bangkok_acc_sp <- as_Spatial(peak_bangkok_acc)
nonpeak_bangkok_acc_sp <- as_Spatial(nonpeak_bangkok_acc)

bangkok_district_sp <- as_Spatial(bangkok_district_format)

peak_bangkok_acc_coords <- coordinates(peak_bangkok_acc_sp)
nonpeak_bangkok_acc_coords <- coordinates(nonpeak_bangkok_acc_sp)

bangkok_bbox_values <- bbox(bangkok_district_sp)

bangkok_window <- owin(xrange = c(bangkok_bbox_values[1, 1], bangkok_bbox_values[1, 2]), 
                   yrange = c(bangkok_bbox_values[2, 1], bangkok_bbox_values[2, 2]))

peak_bangkok_acc_ppp <- ppp(x = peak_bangkok_acc_coords[, 1], y = peak_bangkok_acc_coords[, 2], window = bmr_window)

nonpeak_bangkok_acc_ppp <- ppp(x = nonpeak_bangkok_acc_coords[, 1], y = nonpeak_bangkok_acc_coords[, 2], window = bmr_window)

6.3.3 Handling duplicate points

peak_bangkok_acc_jit <- rjitter(peak_bangkok_acc_ppp, retry = TRUE, nsim = 1, drop = TRUE)

nonpeak_bangkok_acc_jit <- rjitter(nonpeak_bangkok_acc_ppp, retry = TRUE, nsim = 1, drop = TRUE)

Check for duplicates after jittering

any(duplicated(peak_bangkok_acc_jit))

[1] FALSE

any(duplicated(nonpeak_bangkok_acc_jit))

[1] FALSE

6.3.4 Combine point events object and owin object

peak_bangkok_acc_final <- peak_bangkok_acc_jit[bangkok_window]
nonpeak_bangkok_acc_final <- nonpeak_bangkok_acc_jit[bangkok_window]

6.3.5 Kernel Density Estimation (KDE) Analysis

sigma_value2 <- 5000  


peak_kde_bangkok_acc <- density(peak_bangkok_acc_final,  
                       sigma = sigma_value2,  
                       edge = TRUE,  
                       kernel = "gaussian")  

nonpeak_kde_bangkok_acc <- density(nonpeak_bangkok_acc_final,  
                       sigma = sigma_value2,  
                       edge = TRUE,  
                       kernel = "gaussian")  

plot(peak_kde_bangkok_acc, main = "KDE of Accident Data(bangkok peak)")
plot(st_geometry(bangkok_district_format), add = TRUE, border = "black", lwd = 2)

plot(nonpeak_kde_bangkok_acc, main = "KDE of Accident Data(bangkok non-peak)")
plot(st_geometry(bangkok_district_format), add = TRUE, border = "black", lwd = 2)

During peak hours, accident hotspots in Bangkok are mainly concentrated in the central and western districts. The deep yellow and red areas indicate higher accident frequencies, especially in major commuting zones near the city center, such as Phra Nakhon, Bang Rak, and Wang Thonglang. The city center shows a strong concentration of accidents, indicating high traffic pressure, while the western areas, such as Thonburi and Bang Khun Thian, also have significant accident density, reflecting the traffic flow connecting these districts to the city center. In contrast, eastern and northern areas show lower accident density during peak hours, suggesting fewer incidents in these regions.

During non-peak hours, the distribution of accidents becomes more dispersed, though central and western Bangkok remain hotspots. The overall density of accidents decreases compared to peak hours. The central districts continue to be key accident hotspots, particularly along major roads leading to commercial and busy areas. In the western part of the city, accident density also decreases, indicating less traffic pressure during off-peak times. However, accidents are more spread out during non-peak hours, especially in the eastern and western regions, reflecting a broader and more dispersed pattern of incidents.

6.4 K function

From the above analysis, we can identify the high-density districts in Bangkok. We will then select these districts to conduct a K-function analysis, which will help us examine the spatial distribution of accidents within these areas.

6.4.1 import acc data

saveRDS(peak_bangkok_acc,"data/aspatical/peak_bangkok_acc.rds")

peak_bangkok_acc<-readRDS('data/aspatical/peak_bangkok_acc.rds')

6.4.2 intersection data

peak_bangkok_district_acc <- st_intersection(peak_bangkok_acc,bangkok_district_format)

6.4.3 Saphan Sung, Khan Na Yao, Bueng Kum, Bang Kapi, Suan Luang

kfunc_acc_data <- peak_bangkok_district_acc %>%
  filter(ADM2_EN %in% c("Saphan Sung", "Khan Na Yao", "Bueng Kum", "Bang Kapi", "Suan Luang"))

kfunc_boundary <- bangkok_district_format %>%
  filter(ADM2_EN %in% c("Saphan Sung", "Khan Na Yao", "Bueng Kum", "Bang Kapi", "Suan Luang"))

kfunc_road <- st_join(main_bangkok_roads, kfunc_boundary)

lixels_kfunc <- lixelize_lines(kfunc_road,
                         10000,
                         mindist = 5000)         

samples_kfunc <- lines_center(lixels_kfunc)

kfun <- kfunctions(kfunc_road,
                             kfunc_acc_data,
                             start = 0,
                   end = 10000, 
                   step = 1000, 
                   width = 500, 
                   nsim = 50, 
                   resolution = 50,
                   verbose = FALSE, 
                   conf_int = 0.05,
                   agg = 100)

kfun$plotk

From the above plot, in the 0-2750 meter distance, the blue line is significantly above the shaded confidence envelope. This indicates that within this distance, car accidents are spatially clustered in the five selected zones. After this distance, the blue line tends to level off and remains within the shaded area, indicating that beyond 2750 meters, the distribution of accident points becomes more random or follows a pattern similar to spatial randomness. The flattening of the blue line suggests that the clustering effect diminishes as the distance increases, and the accident distribution becomes more spread out and less concentrated.

6.5 Time NKDE for bangkok

6.5.1 acc data filter

peak_bangkok_acc <- peak_acc_utm %>%
  filter(province_en %in% "Bangkok")

nonpeak_bangkok_acc <- non_peak_acc_utm %>%
  filter(province_en %in% "Bangkok")

6.5.2 peak density

peak_nkde_result_bangkok <- nkde(
  lines = lixels_bangkok,                      
  events = peak_bangkok_acc,                     
  w = rep(1, nrow(peak_bangkok_acc)),            
  samples = samples_bangkok,                   
  kernel_name = "quartic",                     
  bw = 500,                                    
  div = "bw",                                  
  method = "simple",                          
  grid_shape = c(200, 200),                    
  verbose = TRUE                               
)

# saveRDS(peak_nkde_result_bangkok,'data/aspatical/peak_nkde_result_bangkok.rds')

6.5.3 non-peak density

nonpeak_nkde_result_bangkok <- nkde(
  lines = lixels_bangkok,                      
  events = nonpeak_bangkok_acc,                     
  w = rep(1, nrow(nonpeak_bangkok_acc)),            
  samples = samples_bangkok,                   
  kernel_name = "quartic",                     
  bw = 500,                                    
  div = "bw",                                  
  method = "simple",                          
  grid_shape = c(200, 200),                    
  verbose = TRUE                               
)

# saveRDS(peak_nkde_result_bangkok,'data/aspatical/peak_nkde_result_bangkok.rds')

# saveRDS(nonpeak_nkde_result_bangkok,'data/aspatical/nonpeak_nkde_result_bangkok.rds')

6.5.4 import data

peak_nkde_result_bangkok <- readRDS("data/aspatical/peak_nkde_result_bangkok.rds")

head(peak_nkde_result_bangkok, 10)

 [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
 [6] 0.000000e+00 3.205377e-06 0.000000e+00 2.432477e-06 0.000000e+00

nonpeak_nkde_result_bangkok <- readRDS("data/aspatical/nonpeak_nkde_result_bangkok.rds")

head(nonpeak_nkde_result_bangkok, 10)

 [1] 0.000000e+00 0.000000e+00 0.000000e+00 2.435628e-06 0.000000e+00
 [6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

6.5.5 visualize data

samples_bangkok$peak_density <- peak_nkde_result_bangkok
lixels_bangkok$peak_density <- peak_nkde_result_bangkok

samples_bangkok$peak_density <- samples_bangkok$peak_density * 10000
lixels_bangkok$peak_density <- lixels_bangkok$peak_density * 10000

samples_bangkok$nonpeak_density <- nonpeak_nkde_result_bangkok
lixels_bangkok$nonpeak_density <- nonpeak_nkde_result_bangkok

samples_bangkok$nonpeak_density <- samples_bangkok$nonpeak_density * 10000
lixels_bangkok$nonpeak_density <- lixels_bangkok$nonpeak_density * 10000

tmap_mode('view')

tm_shape(lixels_bangkok) +
  tm_lines(col = "peak_density", palette = "YlOrRd", title.col = "Density", lwd = 2) +
  tm_shape(peak_bangkok_acc) +
  tm_dots(size = 0.1, col = "blue", alpha = 0.5, title = "Accidents")

tmap_mode('plot')

tmap_mode('view')
tm_shape(lixels_bangkok) +
  tm_lines(col = "nonpeak_density", palette = "YlOrRd", title.col = "Density", lwd = 2) +
  tm_shape(nonpeak_bangkok_acc) +
  tm_dots(size = 0.1, col = "blue", alpha = 0.5, title = "Accidents")

tmap_mode('plot')